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We present high-precision quantum Monte Carlo results for the S — 1/2 XY model on a two- 
dimensional square lattice, in the ground state as well as at finite temperature. The energy, the spin 
stiffness, the magnetization, and the susceptibility are calculated and extrapolated to the thermo- 
dynamic limit. For the ground state, we test a variety of finite-size scaling predictions of effective 
Lagrangian theory and find good agreement and consistency between the finite-size corrections for 
different quantities. The low-temperature behavior of the susceptibility and the internal energy is 
also in good agreement with theoretical predictions. 

PACS numbers: 75.10.Jm, 75.40.Cx, 75.40.Mg 
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I. INTRODUCTION 

Studies of effective continuum field theories have re- 
sulted in detailed predictions for the low-energy-physics 
of quantum spin systems in two dimensions .cTu In or- 
der to make use of these predictions for a given model 
Hamiltonian, the ground state parameters appearing in 
the Lagrangian formulation have to be determined. Spin- 
wave theoryO can in some cases give good estimatesjj but 
in general some numerical method has to be employed 
in order to obtain accurate results. Since the theories 
can also predict the finite-size scaling behavior of vari- 
ous physical quantities, numerical results for a series of 
lattice sizes can be used to extract the ground state pa- 
rameters. Such calculations are also important for testing 
theoretical predictions. 

With todays computers, Lanczos and related exact di- 
agonalization methods .can be used for square lattices 
with up to 6 x 6 spins.LTt] This relatively small maximum 
size, and the small number of different lattices available, 
can make the finite-size scaling procedures problematic 
if sub-leading corrections are significant. It is therefore 
important to consider also alternative methods that can 
reach larger lattice sizes. 

Here we discuss quantum Monte Carlo (QMC) results 
for the S = 1/2 XY model, defined by the Hamiltonian 



SfS?} 



(1) 



where Sf and Sf are the x- and y-components of a spin- 
1/2 operator at site i, and denotes a pair of nearest- 
neighbor sites on a square lattkjCj—.Numerical studies of 
this model have a long history£jll3tl Exact finite-lattice 
calculations gave the first indications that the 0(2) sym- 
metry is spontaneously broken at T = in the thefc 
modynamic limit .Ell This was later proved rigorously.c3 



QMC simulations have psulted in a quite precise, value 
of the ordered momento and have also shownllj that 
there is a Kosterliz-Thojuless (KT) transition at a tem- 
perature T KT « 0.343 J.ll3 For T < T KT the system ac- 
quires quasi-long-range order, i.e., power law decay of the 
spin-spin correlation function and a non-zero spin stiff- 
ness constant. In accordance with the Mermin- Wagner 
theorem JiJ true long-range order develops only at T = 0. 

Our motivation for carrying out calculations for the 
XY model to even higher precision is to test predic- 
tions of effective Lagrangian theories. In particular, the 
chiral perturbation calculations by Hasenfratz and Nie- 
dermayer have resulted in detailed predictions for the 
finite-size and finite-temperature corrections of various 
quantities,!] in some cases beyond leading order. The 
finite-size and finite-temperature scaling behavior of the 
0(3) symmetric Hgisenberg model has been the topic of 
numerous studies ,LrE2rLJilj and the agreement with the 
predictions has been confirmed to high precision. The 
predictions for the 0(2) symmetric XY model have not 
yet been tested exhaustively, however. A recent finite- 
size scaling study of exact energies for systems with up 
to 6 x 6 spins has quite convincingly demonstrated agree- 
ment with the leading finite-size behavior Here we con- 
sider also several other physical observables. This allows 
us to carry out a number of independent tests of the con- 
sistency of the scaling predictions. With access to larger 
system sizes, we can also improve considerably on the ac- 
curacy of the extrapolated ground state parameters and 
finite-size corrections. Our data are sufficiently accurate 
for addressing also sub-leading corrections. In addition to 
calculations in the ground state, we have also studied sys- 
tems at finite temperature, on lattices sufficiently large 
to enable extraction of the leading finite-temperature cor- 
rections in the thermodynamic limit. 

We have used a numerically exact finite-temperatupi 
QMC method based on "stochastic series expansion" , La 
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0.562485(4) 
0.552696(4) 
0.550436(4) 
0.549643(4) 
0.549296(4) 
0.549118(4) 
0.549020(4) 



0.2769(1) 
0.2718(1) 
0.2705(2) 
0.2700(3) 
0.2698(4) 
0.2695(3) 
0.2699(4) 



Ml 
0.13282(2) 
0.11885(4) 
0.1126(2) 
0.1087(2) 
0.1065(3) 



TABLE I. QMC results for the ground state energy, the 
spin stiffness, and the squared magnetization per spin. The 
numbers within parentheses indicate the statistical errors of 
the least significant digit of the results. 
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i.e., importance sampling of diagonal matrix elements of 
the power series expansion of exp(—(3H), where = 1/T 
is the inverse temperature. The ground state can be 
obtained by choosing a sufficiently large j3 value. This 
method has previously been used to study the finite-size 
scaling behavior and the ground .state parameters of the 
square-lattice Heisenberg model. IB Recently, a new way 
of sampling the terms of thepGxpansion was proposed (the 
"operator- loop" algorithm) JH3 which enables calculations 
to even higher precision. We have used this improved 
method to study the ground state of the XY model on 
lattices with N = L x L spins with L up to 16. We have 
carried out calculations at finite temperature for systems 
with L up to 64. 

In Sec. II we discuss the physical quantities that we 
have calculated. The numerical ground state data and 
the results for finite-size corrections and extrapolations 
to infinite size are presented in Sec. III. In Sec. IV we 
discuss calculations at T > 0. Sec. V concludes with a 
brief summary. 

II. CALCULATED QUANTITIES 

The QMC algorithmic can be implemented in two dif- 
ferent ways; with the spin quantization either along the 
z or i axis. The internal energy, E — (H)/N, can be 
easily calculated in both cases. In the z-representation, 
the spin stiffness constant p can also be evaluated. It 
is defined as the second derivative of the internal energy 
with respect to a twist <f> under which, for all nearest- 
neighbor spin pairs in cither the x or y lattice direction, 
the spin-spin interaction is modified according to 

S?Sf+SVS] (SfSf + SfS?) cos (0) 

-S?S}> sin {4>) + SfSJ sin (0) . (2) 

In analogy with the supepfluid density of a boson 
system,Eil one can show thatta 

_ d 2 E{<p) _ {WD + (W„ 2 ) 
P d 2 cf> 2Np ' 1 ' 
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FIG. 1. Ground state energy vs the inverse cubed of the 
system size for L = 8, 10, 12, 14, and 16 (points with error 
bars). The curve is a fit to Eq. urn, including also L = 4 and 
L = 6. 

where E(<f) is the internal energy per spin and W x and 
W y are the "winding numbers", i.e., the net spin cur- 
rents across the periodic boundaries in the x and y di- 
rections that characterize configurations in simulations 
carried out in the ^-representation. We also calculate 
the spin susceptibility, given by 

*-£((£*)")■ 

In the z-representation, it is not convenient to calculate 
spin-spin correlations involving the x or y spin compo- 
nents. We therefore also use an algorithm implemented 
in the x representation, and there calculate the squared 
magnetization Mj; 

^ 2 = ^((e^))- (5) 

Since the 0(2) symmetry is not broken on a finite lattice, 
the magnetization m per spin in the thermodynamic limit 
is given by the infinite-size extrapolated value of 

m = V2M?. (6) 

We use periodic boundary condition and include each 
spin pair only once in Eq. ([lj). For technical details on 
the simulation algorithm we refer to previously published 
workM 



III. GROUND STATE RESULTS 

In order to obtain the ground state, we have carried out 
simulations at inverse temperatures as high as (3 = 1024 
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FIG. 2. Spin stiffness vs the inverse cubed of the system 
size. The line is a fit to the cubic form (^). 



for lattices with L = 4 — 16. In the z-representation sim- 
ulations, we have calculated the ground state energy to 
within relative statistical errors lower than 1CP 5 . The 
relative accuracy of the stiffness is on the order of 10~ 3 . 
Simulations in the ^-representation result in statistical 
errors for the energy roughly twice as large as in the 
z-representation (for simulations of comparable length). 
The squared magnetization is the least accurate quantity, 
with statistical errors of roughly 0.3% for the largest sys- 
tem size considered in this case (L = 12). In Table § we 
list the results used in the finite-size analysis presented 
below. The results for the energy-ior L = 4 and 6 are 
in perfect agreement with previousEifcl exact diagonaliza- 
tions. The magnetization for L — 4 also agrees with the 
exact diagonalization result. However, for L — 6 there 
is an w 5a deviation from the result presented in Rcf. |^. 
We do not know the reason for this disagreement. The 
QMC method has previously been shown to give per- 
fect agreement with isotropic Heisenberg results for both 
L = 4 and L = 6jl3 and therefore a failure for (only) the 
XY-model magnetization for L — 6 would be surprising. 

For a model with 0(2) symmetry, chiral perturbation 
theory predicts the size dependence of the energy asu 



Eq(L) — Eq 



£3 

L ? > 



£5_ 

L 5 



(7) 



with no 0(1/ L ) term. The leading-order correction is 
given by 



e3 



-7C/2, 



(8) 



where the constant 7 = 1.437745 and c is the spin- wave 
velocity. In order to obtain a good fit to our L = 4 — 16 
data, the 0(1/ L 5 ) term has to be included. The fit 
then has a chi-squared value per degree of freedom of 
0.7. Using instead an 0(1/ X 4 ) term gives chi-squared 
w 6. We can therefore conclude that our results sup- 
port the predictionQ that e 4 = 0. In Figure [l] we show 
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FIG. 3. The difference between the energy in the j = 1 and 
j = magnetization sectors graphed vs 1/L 4 . The open circle 
is the exact result for L — 6. The points with error bars are 
the QMC results for L = 8, 10 and 12. The curve is a fit to 
Eq. M}. 



the data for L > 8 along with the best fit. The ex- 
trapolation L —> 00 gives E = -0.548824(2). The 
finite-size correction constants are ea = —0.807(2) and 
eg = —1.07(3). Using Eq. (||) we obtain the spin-wave 
velocity c = 1.123(2). These results are in agreement 
with the previous extrapolations using exact diagonaliza- 
tion data for L = 2,3,4,5,6,0 but the statistical errors 
are considerably smaller. The energy is also in perfect 
agreement with a previous Green's functipa Monte Carlo 
calculation, which gave E = 0.54883(1).^ 

We are not aware of any predictions for the size de- 
pendence of the spin stiffness of the XY model. ForJthe 
Heisenberg model, the leading correction is 0(1/ L)u In 
contrast, our XY data can be very well fit with only an 
0(1 /L 3 ) term; 



r3 



p(L)=p+^ + .. 



(9) 



with the infinite-size value p — 0.2696(2) and the cubic 
correction = 0.47(1). This fit is shown in Figure^. 

The stiffness result is in perfect agreement with the 
value obtained using exact diagonalization and finite-size 
scaling of the ground state energy in higher magnetiza- 
tion sectorsJa The energy per spin in the sector with mag- 
netization j = Y^,i Sf should scale asta 



Sn = ^i + 0(/,l/i 6 ) 



2pL 4 



(10) 



We have also carried out some QMC calculations in the 
j = 1 sector. Figure || shows results for L = 8, 10, and 
12, along with the exact result for L = 6. Using also the 
exact L = 4 result , a fit to Eq. (|l0| ) , including an O (1/L 6 ) 
term, gives c 2 /p = 4.70(4). This value agrees with the 
above separate estimates of c and p. The consistency be- 
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FIG. 4. Magnetization vs inverse system size and a fit to 
the scaling form ([f2|). 



tween the results, obtained in two different ways, clearly 
gives very strong support to the Lagrangian theory. 

The susceptibility can be obtained from c and p using 
the standard hydrodynamic relation 



X = p/<- 



(11) 



which with our values of c and p gives x = 0.2138(8). 
For the magnetization, we assume a finite-size scaling 



Ml{L) = Ml+ a j- + ^ + 



(l2_ 

L 2 



(12) 



A fit to this form, including the quadratic term, is 
shown in Figure ^ and gives M% = 0.0956(6) and a\ = 
0.1212(5). According to Eq. (jq), the magnetization is 
thus m = 0.437(2), in good agreement with a previous. 
Green's function Monte Carlo result [to = 0.441 (5)]. til 
For the Heisenberg model, the linear correction fac- 
tor is related to ground state parameters according to 
ai = aM|/(cx), where a = 0.62075B The oi obtained 
here is instead consistent to within a few percent r-svith 
ai = aM 2 /(2cx)- As in the scaling of the energyja the 
leading size correction is hence proportional to the num- 
ber of gapless modes in the symmetry-broken system. 



IV. FINITE-TEMPERATURE RESULTS 

We now discuss calculations aimed at extracting the 
leading finite-temperature corrections to the ground 
state. Previously, extensive simulations were carried out 
in order to study the KT transitioru-and a critical tem- 
perature Tkt ~ 0.343J was found.t3 Here we focus on 
the susceptibility and the internal energy density at lower 
temperatures. 

Figure || shows the susceptibility, defined in Eq. (||), 
calculated for lattices with L = 32 and 64 for tem- 
peratures T/J > 0.05. Within statistical errors there 
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FIG. 5. Spin susceptibility vs temperature for lattices with 
L = 32 (solid circles) and L = 64 (open circles). The dashed 
line is the estimated constant low-temperature value. The 
inset shows results for L — 32 over a wider temperature range. 



are no differences between the data for the two system 
sizes, and hence the results represent the thermodynamic 
limit. Chiral perturbation theory predicts a temperature- 
independent susceptibility for the XY modelia Within er- 
ror bars, our data is temperature independent for T/J < 
0.15, with the constant value estimated at x = 0.2096(2). 
This result is lower (by about five standard deviations, 
or 2%) than the value extracted in Sec III from the ratio 
X = p/c 2 . The spin-wave velocity was obtained from the 
leading finite-size correction to the ground state energy. 
It is possible that this correction is affected by the pres- 
ence of sub-leading corrections that are not completely 
taken into account by only including up to 0(1/ L 5 ) terms 
in Eq. ([?]) . Using the data shown in Fig. || is a more di- 
rect way of extracting the T = susceptibility, and we 
therefore consider it more reliable. An improved estimate 
of c is then c = p/x = 1.134(2). 

For a finite lattice, the gap between the j = ground 
state and the finite-magnetization states implies an ex- 
ponential decay of x to zero as T — > 0. Chiral pertur- 
bation theory predicts that the form of this decay for 
TL 2 p/c 2 <C 1 is given byQ 



1 



TL 2 



exp[-c 2 /(2pTL 2 )}. 



(13) 



In Figure ^| we show QMC data for L = 4 — 16 along with 
this prediction, where we have used c 2 / p = 4.77, corre- 
sponding to our best estimates of c and p. The agreement 
is not perfect, but satisfactory considering that there are 
no adjustable parameters and that there should also be 
corrections to Eq. (|l3|). For L = 4, we have also calcu- 
lated x using exact diagonalization down to much lower 
temperatures. As shown in Figure 0, there is a very good 
agreement with the predicted form over a sizable low- 
temperature range. Note, however, that the asymptotic 
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FIG. 6. The logarithm of the spin susceptibility for differ- 
ent system sizes vs the logarithm of the temperature (solid 
circles connected by dashed lines, L = 4, 6, . . . , 16 from right 
to left). The solid curves are obtained from Eq. jl^ ) with 
c 2 /p = 4.77. 

T — > decay is always purely exponential, without the 
1/T factor in Eq. (|l5|), as only the lowest j = ±1 states 
contribute. 

Chiral perturbation theory predicts that the Jow- 
temperature form of the internal energy is given byu 

E(T) =E(0) + ^It 3 + O(T 5 ), (14) 

where £(3) = 1.20206. In Figure || we show results for 
T/J > 0.05 calculated for L — 64 lattices. We also show 
a comparison between L = 16, 32, and 64, which suggests 
that L = 64 gives the thermodynamic limit within error 
bars. A fit to the L = 64 data with 0(T 3 ) and 0{T 5 ) 
terms, and including also the T = energy extracted in 
Sec. Ill, gives the cubic correction ((3)/irc 2 = 0.284(5). 
This corresponds to c = 1.16(1), which is higher than the 
value 1.134(2) extracted above using the estimates of p 
and x by about 2.5 error bars. Again, this small deviation 
may indicate some influence of higher-order corrections 
in the energy scaling. 



V. SUMMARY 

In summary, we have presented extensive QMC calcu- 
lations for the two-dimensional 5 = 1/2 XY model. We 
have carried out finite-size and finite-temperature scaling 
of several physical quantities. The results are consistent 
with the predictions of effective Lagrangian theory] to 
within 1-2%. 

The best estimates of the ground state parameters re- 
sulting from our calculations are 
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FIG. 7. The logarithm of the exact spin susceptibility 
for a 4 x 4 lattice vs the logarithm of the temperature 
(dashed curve) compared with the predicted form (b_3|) with 
(? I p = 4.77 (solid curve). 

E Q = 0.548824(2), 
p = 0.2696(2), 
to = 0.437(2), 
X = 0.2096(2), 
c= 1.134(2). 

The ground state energy, E , the spin stiffness, p, and 
the square of the magnetization, to 2 , were all calculated 
directly in the ground state for systems of linear size L = 
4—16 and extrapolated to infinite size. The susceptibil- 
ity x perpendicular to the XY spin-plane was calculated 
at finite temperature and extrapolated to T = 0. Only 
the spin- wave velocity c was obtained by a more indirect 
procedure, using the relation c 2 = p/x- The results arain. 
good agreement pSSjith previous exact diagonalizatiorcj'u 
and QMC work,L3 but the precision is considerably im- 
proved over other estimates. Our results for both the en- 
ergy and the magnetization are in remarkably good agree- 
ment with a series expansion calculation.E3 The magne- 
tization is also in [excellent agreement with second order 
spin-wjave theory.E3 The energy obtained in spin-wave 
theoryO is only 0.2% higher than the numerical result 
obtained here. It would be interesting to carry out high- 
order spin- wave calculations also for the other quantities 
discussed here. 
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FIG. 8. The internal energy calculated for L = 64 lattices, 
and a fit to the form (|l4|). The inset shows data for L = 16 
(squares), 32 (open circles), and 64 (solid circles). 
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